\documentclass[10pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage{graphicx}
\usepackage{lmodern}
\usepackage{fourier}
\usepackage{amsfonts,amsmath,amsthm,amssymb}
\graphicspath{{graphics/}}
\setlength\parindent{0pt}
\usepackage{gensymb}
\usepackage[toc,page]{appendix}
\usepackage[left=2cm,right=2cm,top=2cm,bottom=2cm]{geometry}

\begin{document}
\author{Stephen Duffull, Soroush Safaei, Peter Hunter}
\title{Modelling the Michaelis-Menten System Using Bond Graph Technique}
\date{}
\maketitle

\section{Preamble}

A bond graph is a representation of a system that includes both mass and energy transfers. It differs from most pharmacological representations that consider only mass.\\

The purpose of using this technique is to provide a system that is conserved on mass and energy. We measure mass in mole or molar units and in a closed system (not in equilibrium with the environment) require that the sum of all reactants and products remains constant, such that $\Sigma q=Q$, a constant.  In chemical reactions we are also concerned with internal energy ($U$). If we measure the equilibrium energy of the reactants as the $\Delta U^o_{f_{reactants}}$ and the energy of the products as $\Delta U^o_{f_{products}}$ then conservation of energy in an isolated closed system (where the energy is not in equilibrium with the environment) requires $\Delta U^o_{f_{reactants}}+\Delta U^o_{f_{products}}$ to be constant.

\section{The system}

The Michaelis-Menten approximation of enzyme-catalysed reactions can be represented by:

\begin{equation}
  q_1+q_2\underset{k^{-1}}{\stackrel{k^{+1}}{\rightleftharpoons}}q_3\underset{k^{-2}}{\stackrel{k^{+2}}{\rightleftharpoons}}q_2+q_4
  \label{system}
\end{equation}

In this setting, $q_1$ is the reactant, $q_2$ is the enzyme which is conserved during the reaction, $q_3$ is the intermediate complex, and $q_4$ is the product. It is important to remember in this system that $q_2$ is the same species on both sides of this reaction and therefore is the sum of the enzyme as reactant and enzyme as product. Here $q$ is defined in molar units.\\

The second-order rate constant on from $q_1$ and $q_2$ is $k^{+1}$ (units: $mol^{-1}\, l^{-1}\, s^{-1}$) and the first-order rate constant off is $k^{-1}$ (units: $s^{-1}$), similarly the first-order rate constant from $q_3$ is $k^{+2}$ and the rate constant off is $k^{-2}$. Since formation of the product is generally irreversible then $k^{+2}\gg k^{-2}$.\\

A bond graph of this system will yield a schematic of the energy flow. The bond graphs is:

\begin{figure}[ht]
  \begin{center}
  \includegraphics[scale=0.5]{BG-cs.png}
  \caption{The bond graph representation of Michaelis-Menten closed system.}
  \label{fig:BG-cs}
  \end{center}
\end{figure}

In Figure \ref{fig:BG-cs}, the arrows represent the flow of energy. $Re1$ depicts reaction $1$ and $Re2$ reaction $2$, $q$ has units of $mol\, l^{-1}$ and has the same meaning as in the usual reaction diagram, $\nu$ are velocities of the change in substrate ($mol\, l^{-1}\, s^{-1}$), and $\mu$ the energy per mole ($J\, mol^{-1}\, l^{-1}$). A balanced system is balanced on $\nu$ and $\mu$.\\

The two energy points for conservation are $\mu_3$ and $\mu_2$ which correspond to the generation of $q_3$ and $q_2$ respectively. Note that the recycling of $q_2$ is shown here explicitly as this is required in this closed system.\\

The key to conserving both mass and energy is an attention to units. Here, $q$ has units of  $mol\, l^{-1}$ (molar concentration), $\nu$ has units of $mol\, l^{-1}\, s^{-1}$ (this is also velocity as per MM equation), and $\mu$ has units of $J\, mol^{-1}\, l^{-1}$ (energy per molar concentration).\\

To achieve the purpose, for the system to be conserved for mass and energy, the system must be neutral for mass (in this case $mol\, l^{-1}$), that is $\Sigma q$ is constant, and energy ($J\, s^{-1}$), that is $\frac{\mathrm{d}J}{\mathrm{d}s}=\Sigma \mu \times \nu=-A$ and is constant ($A$ is defined as affinity, when negative the reactants have affinity for each other). This is equivalent to Gibbs free energy at equilibrium for an isolated closed system. A negative value for Gibbs free energy means that a reaction will occur readily.\\
 
It is sufficient in the bond graph representation that only the four points $\mu_2$, $\nu_1$, $\mu_3$ and $\nu_4$ are conserved for energy and mass for the system to meet our goal.\\

\section{The mass system}

Given that the derivative of the moles of substance over time is velocity then $\dot q=\nu$, where $\dot q={\mathrm{d}q}/{\mathrm{d}t}$. Therefore:

\begin{align}
  \begin{split}
    \dot q_1&=-\nu_1,
    \\
    \dot q_2&=\nu_2,
    \\
    \dot q_3&=\nu_3,
    \\
    \dot q_4&=\nu_4.
  \end{split}
  \label{mass-system}
\end{align}

\section{The energy system}

Based on Gibbs free energy and at equilibrium, then the following energy per $mol$ relationships are derived:

\begin{align}
  \begin{split}
    \mu_1&=R.T.ln\left(K_1.q_1\right),
    \\
    \mu_2&=R.T.ln\left(K_2.q_2\right),
    \\
    \mu_3&=R.T.ln\left(K_3.q_3\right),
    \\
    \mu_4&=R.T.ln\left(K_4.q_4\right).
  \end{split}
  \label{energy-system}
\end{align}

Here, $R$ is the gas constant ($=8.314$ $J\, mol^{-1}\, \degree K^{-1}$), and $T$ is the temperature in Kelvin ($=298$ $\degree K$). Finally, $K_x$ is an equilibrium constant and is defined as:

\begin{equation}
  K_x=\left(\dfrac{1}{q_x}\right)e^{\left(\dfrac{\mu_x}{R.T}\right)},
  \label{k1}
\end{equation}
 
and therefore:
 
\begin{equation}
  e^{\left(\dfrac{\mu_x}{R.T}\right)}=K_xq_x,
  \label{k2}
\end{equation}
 
The derivation for Equation \ref{energy-system} is then: 
 
\begin{equation}
  \mu_1=R.T.ln\left(K_1.q_1\right)=R.T.ln\left(e^{\left(\dfrac{\mu_1^0}{R.T}\right)}q_1\right).
  \label{mu}
\end{equation}
 
\section{Conservation of mass and energy}

The velocities are conserved (and therefore mass is conserved) by the two species of interest relating to $\nu_2$ and $\nu_3$. In this model we consider the inward flow of energy and outward flow of energy.  For example, we see that the inward velocity to $\nu_2$ is $\nu_4$ and outward velocity is $\nu_1$, therefore:

\begin{align}
  \begin{split}
    \nu_2&=\nu_4-\nu_1,
    \\
    \nu_3&=\nu_1-\nu_4.
  \end{split}
  \label{cons-mass}
\end{align}

Because the enzyme is reused in the reaction process then the bond graph represents a cycle and hence $\nu_2=-\nu_3$.\\

The energy of the reactants going into the reaction is given by $\mu_1$ (associated with the reactant) and $\mu_2$ (associated with the enzyme). Hence the sum $\mu_5$ is equivalent to $\Delta U_{f_{reactants}}^o$. Similarly the reaction energy of the products $\Delta U_{f_{products}}^o$ is expressed by $\mu_6$ as per Equation \ref{cons-energy}.

\begin{align}
  \begin{split}
    \mu_5&=\mu_1+\mu_2,
    \\
    \mu_6&=\mu_2+\mu_4.
  \end{split}
  \label{cons-energy}
\end{align}

At this point we need to solve for $\mu_1$, $\mu_2$, $\mu_3$, $\mu_4$, $\nu_1$ and $\nu_4$ in order to assess the balance.

\section{Dependence on molar concentration $q$}

It turns out that both velocity and energy are related to molar concentration.  We introduce $\kappa$ as the reaction velocity with units of $mol\, l^{-1}\, s^{-1}$ and recalling that $K$ is our equilibrium constant derived under mass law with units of ${(mol\, l^{-1})}^{-1}$. We see that $K$ is actually the inverse of affinity and is equivalent to ${(K_{eq}^S)}^{-1}$, which is the inverse of the equilibrium binding constant for substrate $S$ in Langmuir’s isotherm.\\

From Equations \ref{mass-system} and \ref{energy-system} and substituting $K$ we can see that:

\begin{align}
  \label{cons1}
  \nu_1=\dot q_1=\kappa_1\left(e^{\left(\dfrac{\mu_5}{R.T}\right)}-e^{\left(\dfrac{\mu_3}{R.T}\right)}\right)=\kappa_1\left(K_1q_1K_2q_2-K_3q_3\right),\\
  \label{cons2}
  \nu_4=\dot q_4=\kappa_2\left(e^{\left(\dfrac{\mu_3}{R.T}\right)}-e^{\left(\dfrac{\mu_6}{R.T}\right)}\right)=\kappa_2\left(K_3q_3-K_2q_2K_4q_4\right).
\end{align}

The link between the reaction rate constants and the thermodynamic parameters are provided by:

\begin{align}
  \label{k+1}
  k^{+1}&=\left(\dfrac{\kappa_1}{q_1+q_2}\right)e^{\left(\dfrac{\mu_1+\mu_2}{R.T}\right)}=\kappa_1 K_1 K_3,\\
  \label{k-1}
  k^{-1}&=\left(\dfrac{\kappa_1}{q_3}\right)e^{\left(\dfrac{\mu_3}{R.T}\right)}=\kappa_1 K_3,\\
  \label{k+2}
  k^{+2}&=\left(\dfrac{\kappa_2}{q_3}\right)e^{\left(\dfrac{\mu_3}{R.T}\right)}=\kappa_2 K_3,\\
  \label{k-2}
  k^{-2}&=\left(\dfrac{\kappa_2}{q_2+q_4}\right)e^{\left(\dfrac{\mu_2+\mu_4}{R.T}\right)}=\kappa_2 K_2 K_4.
\end{align}

We see here that there is a special case solution where the rate constants are specifically drive by the energy transfer since $\mu_3=f(k^{+2},k^{-1})$, here  $f(\, )$ is a function.

\section{The steady state Michaelis-Menten model}

If we assume steady state, where $\nu_2=0$ (the change in abundance of enzyme is constant over time) and $\nu_3=0$ (the change in abundance of intermediate over time is zero), and we also consider that the return rate constant $k^{-2}=0$ then we can define the standard Michaelis-Menten parameters;
the equilibrium constant:

\begin{equation}
  K_m=\dfrac{k^{-1}+k^{+2}}{k^{+1}+k^{-2}}=\dfrac{(\kappa_1+\kappa_2)K_3}{\kappa_1K_1K_2}; \quad \mathrm{units:}\, mol\, l^{-1}=\dfrac{mol\, l^{-1}\, s^{-1}.{(mol\, l^{-1})}^{-1}}{mol\, l^{-1}\, s^{-1}.{(mol\, l^{-1})}^{-2}}
  \label{eqlb_cnst}
\end{equation}

And the maximum velocity:

\begin{equation}
  V_{max}=\kappa_2 K_3(q_2+q_3); \quad \mathrm{units:}\, mol\, l^{-1}\, s^{-1}=mol\, l^{-1}\, s^{-1}.{(mol\, l^{-1})}^{-1}.mol\, l^{-1}
  \label{v-max}
\end{equation}

Here the total abundance of enzyme is $q_2+q_3$. For consistency the units are also solved to ensure mass and energy balance are retained.

\section{The system parameters and laws}

Constitutive parameters: [$\kappa$, $K$]\\
State variables: [$q$]\\
Conditional state variables*: [$\mu$, $\nu$]\\
Mass conservation laws: $\nu_2=\nu_4-\nu_1$, $\nu_3=\nu_1-\nu_4$\\
Energy conservation laws: $\mu_5=\mu_1+\mu_2$, $\mu_6=\mu_2+\mu_4$\\
Prior information: [$k^{-1}$, $k^{+1}$, $k^{-2}$, $k^{+2}$]\\
Constants: [$R$, $T$]\\

*Conditional state variables can be calculated from the state variables, \emph{i.e.} $\nu$ is the time derivative of $q$ and $\mu$ is proportional to $ln(q)$.


\section{The open system}

In order to create an open system from a closed system, the \textit{chemostat} concept is introduced. When one species is fixed to have a constant concentration, a substrate with a fixed state is used as a chemostat. This applies an appropriate external flow to balance the internal flow of the species.\\

Assuming in Equation \ref{system}, $q_1$ and $q_4$ are used as chemostats then we have an open system for the enzyme catalysed reaction. Each chemostat imposes a chemical potential upon the system and inherits the flow on the corresponding junction.

\begin{figure}[ht]
  \begin{center}
  \includegraphics[scale=0.5]{BG-os.png}
  \caption{The bond graph representation of Michaelis-Menten open system.}
  \label{fig:BG-os}
  \end{center}
\end{figure}

In Figure \ref{fig:BG-os}, the dashed box represents a closed system with respect to $q_2$ and $q_3$. However, $q_1$ and $q_4$ correspond to large external pools with effectively fixed concentrations. The only modification that should be made to the closed system to represent an open system is:

\begin{align}
  \begin{split}
    \dot q_1&=0,
    \\
    \dot q_2&=\nu_2,
    \\
    \dot q_3&=\nu_3,
    \\
    \dot q_4&=0.
  \end{split}
  \label{mass-system-os}
\end{align}

The rest of the equations should remain as before.

\subsection{Enzyme degradation}

In this section, an enzyme degradation reaction $Re0$ is added and $\mu_0$ is the zero-potential source. 

\begin{figure}[ht]
  \begin{center}
  \includegraphics[scale=0.5]{BG-dg.png}
  \caption{The bond graph representation of Michaelis-Menten open system with enzyme degradation.}
  \label{fig:BG-dg}
  \end{center}
\end{figure}

\begin{appendices}

\section{OpenCOR output}

The results from OpenCOR for both closed and open system are shown below:

\begin{figure}[ht]
  \begin{center}
  \includegraphics[scale=0.12]{results-cs.png}
  \includegraphics[scale=0.12]{results-os.png}
  \includegraphics[scale=0.12]{results-enzyme.png}
  \caption{From top to bottom: $q_1$ \& $q_2$ \textit{vs} time, $q_3$ \& $q_4$ \textit{vs} time, $v_2$ \& $v_3$ \textit{vs} time, $\mu_5$ \& $\mu_6$ \textit{vs} time. left: closed system, right: open system.}
  \label{fig:results}
  \end{center}
\end{figure}

\section{OpenCOR file}

The CellML file from OpenCOR for both closed and open system is shown below:

\begin{verbatim}
def model MichaelisMenten as
    def comp state as
        var t: second;
        
        // State variables

        var q1: mole {init: 1.0e-0};
        var q2: mole {init: 1.0e-0};
        var q3: mole {init: 1.0e-6};
        var q4: mole {init: 1.0e-6};
        // var v0: mol_per_s; // activate for enzyme degradation
        var v1: mol_per_s;
        var v2: mol_per_s;
        var v3: mol_per_s;
        var v4: mol_per_s;
        var u1: J_per_mol;
        var u2: J_per_mol;
        var u3: J_per_mol;
        var u4: J_per_mol;
        var u5: J_per_mol;
        var u6: J_per_mol;
        // var u0: J_per_mol {init: 0.0}; // activate for enzyme degradation

        // Constitutive parameters
        
        var R: J_per_mol_K {init: 8.324};
        var T: kelvin {init: 300};
        var K_q1: per_mol {init: 20.0};
        var K_q2: per_mol {init: 20.0};
        var K_q3: per_mol {init: 20.0};
        var K_q4: per_mol {init: 20.0};
        var K_Re1: mol_per_s {init: 0.1};
        var K_Re2: mol_per_s {init: 0.1};
        // var K_Re0: mol_per_s {init: 0.1}; // activate for enzyme degradation

        // Conservation laws
        
        ode(q1, t) = -v1; // set to 0 for chemostat
        ode(q2, t) = v2;
        ode(q3, t) = v3;
        ode(q4, t) = v4;  // set to 0 for chemostat
        v2 = v4-v1;       // add +v0 for enzyme degradation
        v3 = v1-v4;
        u5 = u1+u2;
        u6 = u2+u4;
        
        // Constitutive relations
        
        u1 = R*T*ln(K_q1*q1);
        u2 = R*T*ln(K_q2*q2);
        u3 = R*T*ln(K_q3*q3);
        u4 = R*T*ln(K_q4*q4);
        v1 = K_Re1*(K_q1*q1*K_q2*q2-K_q3*q3);
        v4 = K_Re2*(K_q3*q3-K_q2*q2*K_q4*q4);
        // v0 = K_Re0*(exp(u0/(R*T))-K_q2*q2); // activate for enzyme degradation
    enddef;

    def unit per_mol as
        unit mole {expo: -1};
    enddef;

    def unit mol_per_s as
        unit mole;
        unit second {expo: -1};
    enddef;

    def unit J_per_mol as
        unit joule;
        unit mole {expo: -1};
    enddef;
    
    def unit J_per_mol_K as
        unit joule;
        unit mole {expo: -1};
        unit kelvin {expo: -1};
    enddef;
enddef;
\end{verbatim}
\end{appendices}
\end{document}
